
alpha = -180:5:180;
beta = -90:5:90;

alpha_r = deg2rad(alpha);
beta_r = deg2rad(beta);

[alpha_mesh, beta_mesh] = meshgrid(alpha_r, beta_r);

for i1 = 1:size(alpha_mesh,1);
    for i2 = 1:size(beta_mesh,2);
    %R = abs(cos(alpha_mesh).*1);
    alpha_f = alpha_mesh(i1,i2);
    beta_f = beta_mesh(i1,i2);

    alpha_fr = deg2rad(alpha_f);
    beta_fr = deg2rad(beta_f);
    
    beta_fr1 = deg2rad(30);

    R = [ 0.25 -0.25 0 ; -0.25 -0.25 0; -0.25 0.25 0; 0.25 0.25 0];

    K = [ 0 cos(beta_fr1) sin(beta_fr1)];
    K1 = [ cos(alpha_fr)*sin(beta_fr) sin(alpha_fr)*sin(beta_fr) cos(beta_fr)];
    dPhi = 2*pi*K*R';
    B = exp(-1i*dPhi);
    B1 = exp(-1i*2*pi*K1*R');
    R1(i1,i2) = abs(sum(B.*B1));
    end
end
[x, y, z] = sph2cart(beta_mesh, alpha_mesh, R1);

figure(1);
surf( x, y, z)
xlabel('x'); ylabel('y'); zlabel('z')

minc = min( [min(min(x)) min(min(y)) min(min(z))] );
maxc = max( [max(max(x)) max(max(y)) max(max(z))] );
xlim([minc maxc]); ylim([minc maxc]); zlim([minc maxc]);

% surf( alpha_mesh, beta_mesh, R)